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Ai>sfract-The  problem  of  estimating  the  average  Conduction 
Velocity  (CV)  in  the  muscle  fibers  is  considered.  The  velocity  is 
estimated  from  an  array  of  noisy  deterministic  and  unknown 
Motor  Unit  Action  Potentials  (MUAP)  acquired  by  an  array  of  K 
co-linear  electrodes.  In  this  array  a  vector  of  (K-l)  independent 
time  delays  (TD)  are  to  be  estimated.  For  a  given  length  of  the 
array,  the  number  K  has  to  be  determined  in  order  to  improve 
the  estimation  of  the  velocity.  In  this  communication  we  show 
the  influence  of  the  value  of  K  in  the  CV  estimation  in  presence 
of  noise  and  changes  in  the  MUAP  shape  along  the  array. 
Keywords  -  Conduction  Velocity,  array  of  sensors,  change  in 
shape,  Time  Delay  estimation 

I.  Introduction 

The  problem  of  time  delay  (TD)  estimation  is  of  interest  in 
many  applications  such  as  in  Radar,  Sonar  and  in  biomedical 
instrumentation.  Many  authors  tackled  the  problem  of  passive 
TD  estimation  between  two  spatially  separated  sensors  in  the 
presence  of  independent  white  Gaussian  noise.  The  optimal 
method  to  estimate  the  TD  r  is  the  Generalized  Cross- 
Correlation  method  (GCC),  but  this  requires  a  priori 
knowledge  about  signal  and  noise.  In  the  biomedical 
applications,  more  precisely  in  Surface  ElectroMyoGraphy 
(SEMG)  cross-correlation  method  has  been  applied  to 
estimate  the  TD  and  then  the  Conduction  Velocity  (CV)  in 
the  muscle  fibers  [1]:  this  method  could  give  highly  accurate 
estimations  provided  that  the  observed  signals  in  the  array  are 
coherent,  i.e.  neither  scaled  in  time  nor  different  in  shape,  but 
only  delayed  with  an  amount  of  r  between  each  other.  At 
low  Signal  to  Noise  Ratio  (SNR)  or  in  case  of  non  ideal 
acquisition  (e.g.  the  array  is  not  strictly  parallel  to  the  muscle) 
this  condition  is  not  satisfied,  leading  to  high  estimation 
errors.  In  case  of  unknown  signal,  better  TD  estimates  can  be 
obtained  if  the  signal  is  estimated  before  the  correlation 
process.  On  the  other  hand,  when  multiple  realizations  A  of  a 
signal  embedded  in  independent  white  noise  sequences  are 
available,  synchronous  averaging  of  these  realizations 
decreases  the  standard  deviation  of  the  residual  noise  by  an 

order  of  -J~K  .  Consequently,  we  look  forward  to  enhance  the 
TD  estimation  by  applying  synchronous  averaging.  Indeed, 
TD  estimation  by  beamforming  [2],  which  is  itself  a  special 
case  of  the  Maximum  Likelihood  Estimator  (MLE)  of  TD, 
uses  the  synchronous  averaging  to  estimate  the  TD  between 
successive  signals  whose  value  is  assumed  to  be  unique  in  a 
linear  array.  If  the  TDs  between  successive  signals  are  not 
equal,  the  synchronous  averaging  requires  all  the  (A-l) 
independent  TD  estimates  between  the  couples  of  signals  in 
the  array  ;  the  problem  is  to  estimate  all  these  delays. 
Therefore,  if  the  signal  is  unknown  and  the  TDs  are  not  equal, 


the  performances  of  the  TD  estimators  are  decreased.  By 
surveying  the  works  done  in  the  SEMG  domain  and 
especially  the  estimation  of  Conduction  Velocity  (CV)  in  the 
muscle  fibers,  one  can  conclude  that  array  signal  processing 
methods  can  be  used  to  improve  the  CV  estimates.  The  CV  in 
the  muscle  fibers  is  an  important  parameter  of  the  myoelectric 
signal  which  describes  muscle  fatigue  manifestation  during 
voluntary  or  elicited  contractions.  It  may  assess  the 
contraction  level  of  the  muscle,  the  age,  etc..  A  review  of  this 
signal  processing  application,  measurement  techniques  and 
clinical  applications  of  the  CV  are  discussed  in  [3],  In  the 
literature,  many  studies  [1]  have  been  carried  out  using  two 
electrodes  to  acquire  the  necessary  signals  for  the  CV 
estimation.  However,  the  CV  estimation  using  an  array  of 
EMG  signals  was  tackled  by  Schneider  et  al.  [4],  where  the 
estimation  is  carried  out  by  the  inclination  of  the  line  joining 
the  maximums  of  the  MUAPs.  Schneider  showed  that  there  is 
a  remarkable  CV  variability  between  different  sites  along  the 
muscle  fibers  of  a  single  Motor  Unit  (MU),  that  is,  the 
interelectrode  TD  varies  spatially.  This  fact  is  also  proved 
from  theoretical  model  of  SEMG  signals  [5],  therefore,  an 
estimation  of  a  spatial  average  of  the  CV  contains  all  these 
variations.  The  estimated  CV  between  two  successive 
electrodes  v  can  be  found  as  v  =  ll i  where  l  is  the  spatial 
distance  between  the  electrodes  in  meters  and  r  the  estimated 
TD  between  the  two  successive  signals  in  seconds.  The  TD 
estimate  r  is  a  random  variable  (r.v.)  with  an  expected  value 
of  t  .  It  is  obvious  that  uncertainty  on  v  is  a  decreasing 
function  of  the  of  r  value.  Hence,  it  is  better  to  estimate  a 
velocity  from  large  TDs  rather  than  small  ones.  The  question 
is,  why  use  multiple  electrodes  with  small  interelectrode 
distances  while  two  extreme  electrodes  on  the  muscle  unit 
gives  a  lower  variance?  There  are  two  reasons  for  this  : 
firstly,  the  constant  shape  hypothesis  is  more  pertinent 
between  near  signals,  consequently  better  TD  estimation 
using  correlation  ;  secondly,  (it’s  the  aim  of  the  present  work) 
it  can  be  shown  that  the  average  CV  estimation  along  the 
whole  muscle  unit,  can  be  improved  by  using  intermediate 
CV  estimates  under  some  linear  constraints  on  the  TDs.  We 
will  show  in  simulation,  the  influence  of  K  and  changes  in  the 
MUAP  shape  onto  this  improvement. 

II.  Estimation  of  multiple  tds  along  the  array 

Several  models  of  the  recording  of  the  propagating  MUAPS 
along  the  muscle  fiber  has  been  proposed  [6].  We  will  use  the 
following: 

xk{ri)  =  s{n-dk)+wk(ri)  k  =1, ...,  K;n  =  l, ...,  N  (1) 
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where  xk  («)  is  the  signal  recorded  by  the  electrode  number  k. 
In  order  to  simplify  we  can  assume  that  dI  =  0  .  The  wk  («) ’s 
are  mutually  uncorrelated,  with  zero  mean  and  the  same 
variance  a2 .  This  model  is  proposed  assuming  that  the  CV  is 
a  deterministic  variable  that  varies  spatially  (i.e.  along  the 
muscle  unit):  therefore,  the  TDs  between  consecutive  couples 
of  electrodes  will  vary  along  the  array. 

As  shown  in  [6],  the  Maximum  Likelihood  Estimation  of  the 
dk  ‘s  is  given  by: 


1  K- 1  N 

1  A  1 

d  =  arg  min 

- 2x,  (n)  —  Y1xk(n  +  dk-di) 

d 

2(7  *=1  n= 1 

K  k=i 

Defining  =  dt  -  dj ,  the  solution  of  (2)  is  equivalent  to 


„  1  K  K 

T  =  argmax(-— (£fcorr  (tu.)+  E  fcorr  (r2t )  + 

T  (J-K  k= 2  k= 3  (3) 


•  +  fcorr. 


XK-1  >XK  x  L  ( A'-l)  A 

where  fcorr(  )  stands  for  the  cross-correlation  function  and 
defining  vectors  as  : 

d  =  [dl,d2,...,dKf ,  xk=[xk(l),xk(2  ),...,xk(N)f  and 

A  j.  A  A  A  A  A  A  -.'J' 

"T  =  LTi2>  Tij,..  Tjaj  *23’  ■  ■  ■  ’  T2K  ’  ■  ■  ■  ’  K-1)K  1 

Thus  the  optimal  result  is  given  by  performing  an  exhaustive 
search  over  the  d: ’s,  a  (K-\)  dimensional  maximization,  with 


substituting  the  r,( by  dj-ddm  (3).  Indeed  the  method  is 

optimal  but  time  consuming  when  high  resolution  is  needed. 
An  alternative  approach  uses  a  suboptimal  result  given  by  a 
post  correlation  processing:  maximizing  each  term  of  (3)  with 
respect  to  all  shift  operations,  we  obtain  the  correlation 
estimates  Tij  ,we  call  coarse  estimates.  Note  that  the  TDs  are 

estimated  continuously  by  fitting  a  parabola  on  the  apex  of 
the  cross-correlation  [1].  Then  exploiting  the  K(K- 1)/2  linear 
relations  between  the  as  constraints  (Chasles  relations), 

leads  to  fine  estimates.  The  advantage  of  this  approach  is  the 
rapidity  and  high  resolution  result  with  respect  to  the  optimal 
approach.  This  suboptimal  approach  can  be  presented  as 
following. 

First  let  us  define  the  K(K-l)/2  linear  relations  between  the 
tjj .  These  linear  independent  relations  (Chasles  relations) 
will  supply  additional  information  and  can  be  exploited  to 
improve  the  t  estimate.  These  relations  can  be  written  as: 


712  +  T23  —  T13 


T12  +  r2A  T1A 
T23  +  T34  =  T24 


t23  +  Ua  —  t2K 
r(A-2)(A-l)  +  T(A-1)A  =  T(K-2)K 


In  matrix  form  this  can  be  written  as: 

Ct  =  0 

where  C  is  the  matrix  of  constraints  with  elements  (-1,0,1) 
and  size  of  (K-l)(K-l)/2xK(K-l)/2  and  has  a  full  rank.  This 
matrix  of  constraints  is  involved  in  the  final  estimation  step 
consisting  in  minimizing  an  error  vector  e  in  the  set  of  linear 
equations  that  relates  the  vector  to  be  estimated  t  to  the 
continuous  coarse  TD  already  estimated  t  together  with  the 
constraints: 


"o" 

C 

"o 

= 

T+ 

T 

i 

e 

where  0  is  a  ( K-l)(K-2)/2xl  zero  vector,  I  is  a  K(K-\)/2xK(K- 
l)/2  identity  matrix,  t  is  the  K(K- l)/2xl  vector  of  coarse 
interelectrode  TD  estimates.  The  vector  e  contains  K(K-l)/2 
error  elements  with  e-  =  r;j  -  rs  . 


The  least  square  estimator  of  the  TDs  t  from  (4)  and  defining 


,  would  be  given  by: 


t  =  (QR  >Q)  >QTR  1 


0 

T 


(5) 


with  R,  the  covariance  matrix  of 


.  Using  some  particular 


LeJ 

assumption  one  gets  the  final  result: 

T  =  [l-®CT(COCT)  ‘c]t  (6) 

with  <D  =  yl  and  y  the  variance  of  the  error  e.  Note  that  due 


to  the  form  of  (6)  the  value  of  y  is  not  needed. 

Since  xis  given  by  a  suboptimal  approach,  the  variance  of 
the  ij  cannot  reach  the  Cramer-Rao  Lower  Bound  (CRLB) 
given  by  [6],  for  any  value  of  K: 

var(  r,( )  =  varL/  ; -d,  )  =  var (dt )  > 

where  s’ is  the  first  derivative  of  the  signal  with  respect  to 
time.  One  could  notice  that  this  theoretical  bound  doesn’t 
depend  on  the  number  of  sensors  (K)  meaning  that  increasing 
K  doesn’t  lead  to  a  performance  improvement.  In  fact,  the 
estimator  variance  will  tend  toward  the  CRLB  as  K  tends  to 
infinity  meaning  that  the  accuracy  of  the  proposed  estimator 
will  be  a  function  of  K. 

This  improvement  of  the  TDs  estimation  is  crucial  since  the 
CV  and  TDs  are  linked  by  the  inverse  function.  As  we  will 
see  in  the  following,  the  calculation  of  the  CV  using  the  set  of 
TDs  can  be  achieved  in  various  way. 


III.  Velocity  estimation 

The  velocity  can  be  estimated  from  the  well  known  relation 
v  =  /  /  d  where  l  is  the  distance  between  two  consecutive 
electrodes  and  d  is  the  TD  between  these  electrodes  expressed 
in  number  of  samples.  An  analytic  form  of  the  pdf  of  the 
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velocity  estimator  can  be  obtained  provided  that  the  TD  pdf  is 
known.  The  choice  of  Gaussian  pdf  is  a  good  one  and  can  be 
satisfied,  since  the  TD  estimators  are  obtained  from  cross¬ 
correlations  and  this  is  a  MLE.  Moreover,  the  noise  processes 
are  considered  to  be  Gaussian.  The  relation  between  the  delay 
pdf  p-  and  the  velocity  pdf  p  -  is  given  by  : 


Py(v)  =  4ft(“) 

V  V 

where  p  =  lfs,  fs  being  the  sampling  frequency  in  Hertz 
and  1  is  the  distance  in  meters.  Taking  a  Gaussian  pdf  for  the 
fine  TD  estimators  r  ,  with  mean  p-  (this  is  equal  to  the 

exact  TD  for  a  non-biased  delay  estimator)  and  variance  cs\ , 
then  the  velocity  pdf  is  given  by  : 


Py(v)  = 


P 

V^rcrjV2 


exp 


(p-vPif 

2a:v2 


(7) 


One  should  note  that  from  a  theoretical  point  of  view  the 

mathematical  expectation  of  V  using  (7)  as  pdf  doesn’t  exist. 
Using  truncation  in  numerical  calculation  avoid  this  problem. 
For  a  given  velocity  (i.e.  a  given  constant  ratio  of  // /rfwith 
arbitrary  1  and  pt  values)  the  variance  of  the  velocity 
estimator  will  be  high  when  1  and  p.  are  small,  because 
noting  the  expression  of  pAv) ,  we  can  see  that  the 


mathematical  expectation  of  V  explicitly  depends  on  /  and 
p-  and  not  only  on  the  ratio  Up,  -  This  remark  is  illustrated 


2 

by  giving  mean  ( p- )  and  variance  ( <r. )  of  V  for  different 


values  of  on  l  and  pi  with  a  constant  ratio  of  the  two 
variables,  as  in  Tab.l. 


Table i 


MEAN  AND  VARIANCE  OF  CV 


Pi 

/ 

Pv 

2 

G  „ 

V 

5 

0.01 

4.16 

18.3 

10 

0.02 

4.04 

16.5 

15 

0.03 

4.01 

16.2 

From  the  precedent  theoretical  analysis  we  conclude  that  the 
smaller  the  distance  between  electrodes  (smaller  TD),  the 
higher  will  be  the  velocity  estimator  variance.  This  leads  us  to 
estimate  TDs  from  widely  spaced  electrodes,  but  here  again 
we  are  limited  from  the  muscle  unit  length  and  spatial 
velocity  resolution  points  of  view.  Another  reason  to  take  an 
array  of  signals  with  small  interelectrode  distance,  is  to  get 
shape  resemblance  between  signals  in  the  array  and 
consequently  better  TD  estimation  using  correlation. 

Keeping  in  mind  that  we  have  a  vector  of  TD  estimates 
obtained  from  different  couple  of  electrodes  in  the  electrode 
array,  it  is  expected  that  an  average  estimate  of  the  velocity 
from  all  the  independent  TDs  will  decrease  the  velocity 
variance.  The  average  velocity  is  the  harmonic  average  vH , 


since  we  have  equal  electrode  spacing  and  probably  varying 
velocity.  This  is  calculated  as  : 


1 


1 


V 


v 


K- 1 


A 


1  1 

— +  — + 
v,  v. 


1 


where  v;  =  lfs  /  ii_li ,  then 


v„  = 


(K~WS 


^  T12  +  T23  +  •••  +  j 


(8) 


When  the  total  length  of  the  sensors  array  is  given  equal  to  L, 
and  adapted  to  the  length  of  the  muscle  fibers,  the  distance 
between  successive  sensors  is  1  =  L/(K- 1) .  So,  (9) 
becomes: 

f  ,  „  A 


Lfs 


y  12  +  T23  +  '  • '  +  Z{K-l)K  j 


(9) 


The  average  CV  can  be  also  estimated  by  taking  the  TD  (  rL ) 
(between  the  two  most  separated  signals  (using  the  most 
separated  sensors)  only: 


(10) 


U 

But  as  mentioned  above  there  are  some  disadvantages  when 
used  with  real  signals  because  of  the  signal  shape  variation 
between  larger  spaced  electrodes.  From  a  theoretical  point  of 
view  and  taking  the  CRLB  as  variance  for  the  estimated  r tj , 


the  variance  of  the  sum  rl2  +  r23  + . . .  +  r{K_1)K  (appearing  in 


(9))  is  increasing  as  K  becomes  larger,  leading  to  a  loss  of 
accuracy  in  the  CV  estimation.  This  fact  seems  to  be  in 
contradiction  with  the  assumption  that  the  accuracy  of  CV 
estimation  will  be  higher  as  the  number  of  sensors  increases. 
We  will  see  in  the  following  that  since  practically  the 
variance  of  r tj  is  not  the  CRLB  one,  this  conclusion  is  not 
valid. 


IV.  Simulations  and  results 

In  the  following  simulations,  we  will  show  that  on  the 
contrary  to  the  ideal  approach  (TDs  variance  equalizes  the 
CRLB)  the  performance  of  CV  estimation  (9)  is  improved 
when  taking  large  value  of  K.  We  also  show  the  influence  on 
the  estimates  (9)  and  (10)  of  a  change  in  shape  of  the  MUAP 
along  the  fiber  when  K  increases.  In  both  cases  the  duration 
of  the  simulated  signal  is  about  40  samples  (corresponding  to 
a  20ms  duration  MUAP  sampled  at  a  2000Hz  sampling  rate) 
and  the  constant  CV  equal  1  m/sec.  The  criteria  used  in  order 
to  characterize  the  accuracy  is  the  Mean  Square  Error  (MSE) 

defined  as  mse  =  a\  +  bias1  since  it  has  been  shown  that 
uncertainty  in  the  TDs  estimation  adds  bias  in  the  CV 
estimation. 

A.  Influence  ofK  value  on  the  CV  estimation 
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Fig.  1:  MSE  of  CV  estimation  (9)  and  (10)  for  K=3,5,10,20  (circle)  and  K=2 
(only  two  sensors )( stars) 

In  this  simulation  the  SNR  is  about  18dB  and  the  number  of 
trials  to  calculate  the  mse  is  equal  to  500.  The  result  of  this 
simulation  is  given  in  Fig.  1  where  the  mse,  characterizing  (9) 
and  (10),  defined  above  is  a  function  of  K  =  (3,5,10,20) . 

B.  Influence  of  change  in  shape  on  the  CV  estimation 

The  aim  of  this  simulation  is  to  show  the  advantage  in  using 
more  than  two  sensors  when  the  signal  is  subject  to  change  in 
shape  along  the  muscle  fiber.  In  Fig.2,  signals  used  in  this 
simulation  (no  noise)  are  shown  for  K=  3,  intermediate  shape 
are  easily  deduced  when  K  is  higher  than  3. 

In  Fig. 3,  mse  for  (9)  and  (10)  are  calculated  adding  a  change 
in  shape.  Since  signal  are  not  noisy,  mse  criteria  reveals  only 
the  bias. 


Fig.  3:  MSE  of  CV  estimation  for  K=3,5,10,20,30  (circle)  and  K=2  (only  two 
sensors )( stars) 

V.  Discussion 

In  this  communication,  in  order  to  calculate  the  Conduction 
Velocity  in  the  muscle  fiber.  Time  Delays  are  estimated  from 
an  array  of  K  sensors.  After  giving  a  theoretical  study  of  the 
CV  estimation  performances,  we  have  shown  in  simulations 
the  obvious  advantage  in  using  more  than  2  sensors  (K> 2)  but 
due  to  the  physical  limitation  (sensors  size)  and  the  limited 
gain  (Fig.  land  Fig. 3)  for  A>15,  a  bound  for  the  increasing  of 
K  can  be  proposed.  When  focusing  on  the  bias  reduction  in 
Fig. 3,  the  improvement  is  not  very  important  probably  due  to 
the  signal  type  used  for  the  simulation. 
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